library(matrixStats)
library(sp)
library(maptools)
## Checking rgeos availability: TRUE
library(maps)
library(raster)
library(colorRamps)
library(rgdal)
## rgdal: version: 1.2-6, (SVN revision 651)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-4
library(scales)
library(animation)
library(dismo)
library(deldir)
## deldir 0.1-14
library(rgeos)
## rgeos version: 0.3-23, (SVN revision 546)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Linking to sp version: 1.2-4 
##  Polygon checking: TRUE
library(gstat)

Organizing Breeding Timing Data

Here I transform the breeding timing data into a SpatialPointsDataFrame and then create a presence/absence matrix to show which birds are breeding on a given date.

load("Global_bird_breeding_dates_julian.Rdata")

# change dataframe into SpatialPointsDataFrame
spring_data <- data.frame(spring_data)
class(spring_data[,5]) <- "numeric"
class(spring_data[,6]) <- "numeric"
nonas <- spring_data[which(is.na(spring_data[,5])==FALSE | is.na(spring_data[,6])==FALSE),]
nonas <- nonas[which(is.na(nonas[,5])==FALSE & is.na(nonas[,6])==FALSE),]

coord <- coordinates((nonas[,6:5]))
bird_pts <- SpatialPointsDataFrame(data =nonas, coords = coord, proj4string = CRS("+proj=longlat +datum=WGS84"))

# make a presence absence matrix
namevector <- vector(length = 365)

for (i in 1:365){
  namevector[i] <- paste0("day", i)
}

dates_matrix <- data.frame(matrix(ncol = 365, nrow = 1780))
dim(dates_matrix)

colnames(dates_matrix) <- namevector
dates_matrix[] <- 0

bird_pts_binded <- spCbind(bird_pts, dates_matrix)
bird_pts <- bird_pts_binded

i <- namevector[1]
for (i in namevector){
  bird_pts[which(bird_pts$Egg_start_julian_date <= as.numeric(substr(i, 4, nchar(i))) & as.numeric(substr(i, 4, nchar(i))) <= bird_pts$Fledge_end_julian_date), 42 + as.numeric(substr(i, 4, nchar(i)))] <- 1
  bird_pts[which(bird_pts$Egg_start_julian_date > bird_pts$Fledge_end_julian_date), 42 + as.numeric(substr(i, 4, nchar(i)))] <- 3
  bird_pts[which(bird_pts[[i]] == 3 & as.numeric(substr(i, 4, nchar(i))) >= bird_pts$Egg_start_julian_date), 42 + as.numeric(substr(i, 4, nchar(i)))] <- 1
  bird_pts[which(bird_pts[[i]] == 3 & as.numeric(substr(i, 4, nchar(i))) <= bird_pts$Egg_end_julian_date), 42 + as.numeric(substr(i, 4, nchar(i)))] <- 1
  bird_pts[which(bird_pts[[i]] == 3), 42 + as.numeric(substr(i, 4, nchar(i)))] <- 0
}

bird_pts$day1 <- 0
plot(bird_pts[which(bird_pts$day1 == 1),])

save("bird_pts", file = "bird_points_presence_absence.Rdata")
load("bird_points_presence_absence.Rdata")

Plotting Breeding Birds Pre-Climate Change and Post-Climate Change

Here I plot points representing breeding birds throughout the year on top of the 1941-1950 raster and on top of the 2007-2016 raster.

list_dates <- seq(as.Date("2017-01-01"), as.Date("2017-12-31"), by = "day")
list_dates <- format(list_dates, format = "%b %d")

list_files_south <- list.files("Past Interpolation Southern", full.names = TRUE)
list_files_north <- list.files("Past Interpolation Northern", full.names = TRUE)

for (i in 1:365){
  list_files_south[i] <- paste0("Past Interpolation Southern/" , i, "past_northern.Rdata")
  list_files_north[i] <- paste0("Past Interpolation Northern/" , i, "past_northern.Rdata")
}

list_files_south_mod <- list.files("Present Interpolation Southern", full.names = TRUE)
list_files_north_mod <- list.files("Present Interpolation Northern", full.names = TRUE)

for (i in 1:365){
  list_files_south_mod[i] <- paste0("Present Interpolation Southern/" , i, "present_southern.Rdata")
  list_files_north_mod[i] <- paste0("Present Interpolation Northern/" , i, "present_northern.Rdata")
}

c1 <- colorRampPalette(c("navy", "blue", "magenta2", "indianred1", "yellow", "white"), bias = 8.5)

saveGIF({
  for (i in 1:365){
    par(mfrow=c(1,2))
    load(as.character(list_files_south[i]))
    plot(masked, main = paste0("1941-1950 Growing Degree Day Accumulation | ", list_dates[i] ),  mar = c(2, 4, 1, 2) +0.1, col = c1(5000), breaks = seq(0,5000, by = 100), legend = FALSE, xlim = c(-180, 180), ylim = c(-90,90), asp = 1)
    abline(h = 0, lty = 2)
    load(as.character(list_files_north[i]))
    plot(masked, add = TRUE, col = c1(5000), breaks = seq(0,5000, by = 100), legend = FALSE)
    plot(bird_pts[which(bird_pts[[paste0("day", i)]]==1),], add = TRUE, pch = 16, cex = 0.5, col = "lawngreen")
    
    load(as.character(list_files_south_mod[i]))
    plot(masked, main = paste0("2007-2016 Growing Degree Day Accumulation | ", list_dates[i] ),  mar = c(2, 4, 1, 2) +0.1, col = c1(5000), breaks = seq(0,5000, by = 100), legend = FALSE, xlim = c(-180, 180), ylim = c(-90,90), asp = 1)
    abline(h = 0, lty = 2)
    load(as.character(list_files_north_mod[i]))
    plot(masked, add = TRUE, col = c1(5000), breaks = seq(0,5000, by = 100), legend = FALSE)
    plot(bird_pts[which(bird_pts[[paste0("day", i)]]==1),], add = TRUE, pch = 16, cex = 0.5, col = "lawngreen")
    print(i)
  }
}, movie.name = "global_interpolation_birds.gif", ani.width = 1000, ani.height = 500, interval = 0.3)

Extracting Accumulated Degree Days

Here I extract the number of degree days that have accumulated by the day that the birds lay their first eggs. I extract values from the pre-climate change raster and from the post-climate change raster.

1941-1950 Raster Extractions

nonas <- bird_pts[which(is.na(bird_pts$Egg_start_julian_date) == FALSE),]

list_files_south <- list.files("Past Interpolation Southern", full.names = TRUE)
list_files_north <- list.files("Past Interpolation Northern", full.names = TRUE)

for (i in 1:365){
  list_files_south[i] <- paste0("Past Interpolation Southern/" , i, "past_northern.Rdata")
  list_files_north[i] <- paste0("Past Interpolation Northern/" , i, "past_northern.Rdata")
}

# only looking at the northern hemisphere
nh_pts <- nonas[which(nonas$Lat > 0),]

start_dd <- data.frame()

for (i in 1:length(nh_pts)){
  x <- round(nh_pts@data[i,35])
  if(x > 0){
  load(as.character(list_files_north[x]))
  start_dd[i,1] <- extract(masked, nh_pts[i,])
  }
  else{
    start_dd[i,1] <- NA
  }
}

# median degree days accumulated on day of first egg
median(start_dd[,1], na.rm = TRUE)
## [1] 26.6932
# mean degree days accumulated on day of first egg
mean(start_dd[,1], na.rm = TRUE)
## [1] 315.9695

2007-2016 Raster Extractions

list_files_south_mod <- list.files("Present Interpolation Southern", full.names = TRUE)
list_files_north_mod <- list.files("Present Interpolation Northern", full.names = TRUE)

for (i in 1:365){
  list_files_south_mod[i] <- paste0("Present Interpolation Southern/" , i, "present_southern.Rdata")
  list_files_north_mod[i] <- paste0("Present Interpolation Northern/" , i, "present_northern.Rdata")
}

# only looking at the northern hemisphere
nh_pts <- nonas[which(nonas$Lat > 0),]

start_dd <- data.frame()

for (i in 1:length(nh_pts)){
  x <- round(nh_pts@data[i,35])
  if(x > 0){
  load(as.character(list_files_north_mod[x]))
  start_dd[i,1] <- extract(masked, nh_pts[i,])
  }
  else{
    start_dd[i,1] <- NA
  }
}

# median degree days accumulated on day of first egg
median(start_dd[,1], na.rm = TRUE)
## [1] 136.6719
# mean degree days accumulated on day of first egg
mean(start_dd[,1], na.rm = TRUE)
## [1] 655.3654